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We study theoretically and numerically how hard frictionless particles in random packings can 
rearrange. We demonstrate the existence of two distinct unstable non-linear modes of rearrange- 
ment, both associated with the opening and the closing of contacts. Mode one, whose density is 
characterized by some exponent 9' , corresponds to motions of particles extending throughout the 
entire system. Mode two, whose density is characterized by an exponent 9' , corresponds to the 
local buckling of a few particles. Mode one is shown to yield at a much higher rate than mode two 
when a stress is applied. We show that the distribution of contact forces follows P{f ) ~ .*)^ 
and that imposing that the packing cannot be densified further leads to the bounds 7 > j^!^/ and 
7 > where 7 characterizes the singularity of the pair distribution function g{r) at contact. 

These results extend the theoretical analysis of [l] where the existence of mode two was not consid- 
ered. We perform numerics that support that these bounds are saturated with 7 « 0.38, 6 ~ 0.17 
and 9' « 0.44. We measure systematically the stability of all such modes in packings, and confirm 
their marginal stability. The principle of marginal stability thus allows to make clearcut predictions 
on the ensemble of configurations visited in these out-of-equilibrium systems, and on the contact 
forces and pair distribution functions. It also reveals the excitations that need to be included in a 
description of plasticity or flow near jamming, and suggests a new path to study two-level systems 
and soft spots in simple amorphous solids of repulsive particles. 



I. INTRODUCTION 

The dynamics in amorphous materials is so slow that 
thermal equilibrium cannot be reached. In these sys- 
tems properties are history- dependent, and configura- 
tions of equal energy are not equiprobable. Understand- 
ing the properties of these materials is thus a challenge, 
as it requires some knowledge on the dynamics and the 
configuration space it explores. It appears that one 
key aspect of amorphous solids is the possibility for a 
group of particles to rearrange locally. This phenomenon 
is believed to govern the low-temperature properties of 
glasses, where the excitations governing the specific heat 
are not phonons but two-level systems where a group 
of particles can switch between two distinct configura- 
tions [2]. Local rearrangements, the so-called shear- 
transformation zones, also dominate the rheological prop- 
erties of amorphous solids under shear [3H5]. Near the 
glass transition where structural relaxation is thermally 
driven, is it also observed that relaxation occurs in fa- 
vored locations [H [7]. These soft spots correspond to 
regions where low-frequency vibrational modes are quasi- 
localized [S^. In all these situations elucidating the na- 
ture of these local rearrangements, and determining what 
controls their density, has remained a challenge. 

Here we study these questions in packings of hard fric- 
tionless spheres, perhaps the simplest model of amor- 
phous solids. It is widely used as a model system for 
structural glasses and for granular systems, and its micro- 
scopic structure has been studied carefully [9Hl3] . It is 
found in particular that in such packings there are many 
particles that are almost touching, but not quite: the dis- 
tribution function of the gaps between particles g{h) has 
a singularity near contact g{h) ~ h~'^ , with 7 « 0.4 [TOl - 
[T^ . It has been recently noticed that the distribution 



of the contact force amplitude P{f) is also character- 
ized by a non trivial exponent P{f) ^ [T^J US] where 
9 e [0.18,0.45], depending on the system preparation. 
Simple toys model [HI [TBi that have been proposed early 
on to compute P{f) give 9 = 0. In general, comput- 
ing the structure is complicated because it requires some 
knowledge on the ensemble of configurations visited by 
the dynamics. Propositions considering that all mechani- 
cally stable states are equally likely [TB] , or that packings 
are maximally random |17) have so far not been able to 
rationalize these findings. Mean- field approaches based 
on the replica method |12l HH] make quite accurate esti- 
mations on the range of values at which jammed packings 
can be found, but currently predict = 7 — 0. 

Here we shall contend that both the structural prop- 
erties of packings and their low-energy excitations can 
be understood in a common framework, based on the 
principle of marginal stability. This principle, illus- 
trated in Fig.[T] states that the jammed configurations at 
which the system arrives are just stable toward the relax- 
ation processes that drive the dynamics in the unjammed 
phase. If realized, this principle is powerful, both because 
it reduces the configurational space of the jammed states 
encountered, and because it implies an abundance of low- 
energy excitations that are expected to strongly affect 
the response of the system. The principle of marginal 
stability has been successfully applied in some glassy 
systems with long-range interactions, most importantly 
in Coulomb glasses where it implies that the density of 
states must vanish at the Fermi energy [I9l - f23] . but also 
in fully-connected spin glasses P^H^ where it states that 
the distribution of local fields must vanish at vanishing 
local fields [24,. In both cases, marginality leads to the 
presence of power-law avalanches of rearrangements un- 
der forcing, referred to as crackling noise |27j . 
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FIG. 1: Illustration of a marginal stability diagram. A and B 
are observables characterizing the structure. The thick black 
line corresponds to marginal stability: it separates a region 
where dynamical modes or excitations are stable, from a re- 
gion where these modes are unstable. The arrow-decorated 
blue line illustrates a dynamical trajectory occurring when 
the system is prepared. The system is initially unstable, until 
it reaches the marginality line. At this point if the tempera- 
ture or the tapping that drives the system is sufficiently small, 
the dynamics will stop, and the system will lie close to the 
marginality line. In the case of soft spheres, modes are vi- 
brational, A corresponds to the coordination z and B to the 
compressive strain e ~ </> — </)c [281 129| , where <j!> is the packing 
fraction. In the present work modes are non-linear excitations 
associated with opening and closing contacts, A corresponds 
to the exponent 6 or 8' characterizing some aspects of the 
force distribution, and B is the exponent 7 describing the 
distribution of gaps between particles that are not in contact. 
On the stability line many soft excitations are present, and 
rich dynamics, such as crackling noise, can occur. 



Marginal stability has been applied previously to com- 
pressed packings of soft particles [2B1 !M] and to colloidal 
glasses [301 [SB- In these cases linear stability is con- 
sidered, and the modes at play are simply vibrational. 
Requiring that packings lie close to an elastic instability 
yields z — Zc ~ — (jjc]^^'^ where z is the coordination, Zc 
is twice the spatial dimension, and |(/)— (/)c| is the distance 
from the jamming threshold, expressed in terms of pack- 
ing fraction. This scaling is indeed observed numerically 
[91 |3U] and in emulsions [32] . It can be shown using vari- 
ational arguments [5^1 [33] or effective medium approx- 
imation [HU [35] that this behavior of the coordination 
implies various elastic anomalies and two associated di- 
verging length scales near the jamming threshold, which 
are indeed observed [35HiD] . However, this approach is 
mostly limited to linear properties and cannot explain 
some aspects of the structure, e.g. the distribution of 
forces, and leads to a limited insight on the non-linear 
mechanisms that govern plastic flow [HJ [35] . 

To bridge this gap, one of us has considered relaxation 
mechanisms associated with changes of contact networks 
in packings [I] . The excitations considered were extended 
in space, and imposing their stability led to an inequality 
between the exponent characterizing the distribution of 
contact forces and the exponent characterizing the dis- 



tribution of gaps between particles not in contact. It was 
assumed that the density of these excitations is controlled 
by the distribution of contact forces at low-force. Here 
we revisit this approach and incorporate significant mod- 
ifications. Most importantly, we find that they arc two 
distinct unstable excitations associated with the open- 
ing and closing of contacts, instead of one. The novel 
excitations correspond to local rearrangements of a few 
particles. This finding is important, both because local 
excitations are expected to affect various properties of 
amorphous solids as discussed above, but also because, 
at least for the system preparation we use, these local ex- 
citations appear to be more numerous that the extended 
excitations introduced in [1], and therefore govern the 
distribution of contact forces. 

This paper is organized as follows. In Section |ll] the 
relationship between weak contact forces and packing ge- 
ometry is investigated, and in Section jlllj the nonlinear 
stability criterion of such contacts toward opening is de- 



rived. In Section IV a scaling analysis of this criterion 
is performed, which reveals the presence of two relax- 
ation mechanisms, and yields bounds on their respective 
densities. In Section [V] we test numerically our theoreti- 
cal predictions and establish the existence of two distinct 
modes of relaxation. We also show that jammed pack- 
ings produced under our protocol are indeed marginally 
stable. We demonstrate how our findings imply relations 
between P{f) and g{h). In Section VI we argue that 



marginal stability leads to a simple explanation for the 
power-law avalanches of plasticity observed in packings 
under loading [H], and in Section VII we discuss how 
our arguments extend to soft compressed particles. Sec- 
tion IVIIII summarizes our results and indicates relevant 
open questions. 



II. GEOMETRIC NATURE OF SMALL FORCES 
IN PACKINGS 

Contacts carrying small forces play a special role in 
packings of hard particles, as they are the most likely to 
open and lead to rearrangements when external stresses 
are applied. The density of weak forces thus strongly 
affects how the solid can how, and it is important to 
understand what causes certain contacts to carry weak 
forces. Contact forces can be expressed geometrically 
in random packings of frictionless particles, as we now 
recall. We shall initially assume that the packing is con- 
tained in a cubic box of volume made of rigid walls, and 
present the main theoretical results in this context. Our 
results are straightforward to extend to periodic bound- 
aries, that we shall use to perform numerical tests. We 
consider a disordered packing of N hard frictionless par- 
ticles of diameter uq, in spatial dimension d. The pack- 
ing is formed by pushing particles together by reducing 
the box size, so as to apply a pressure p that fixes the 
scale of contact forces in the system. Microscopically, the 
boundaries apply external forces Fi on all the particles i 
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in contact with it. We denote the i"" particle's position 
by Ri, the pairwise vector of differences by Rij = Rj — Ri, 

and the bare pairwise distance by = \J R^j ■ Rij . 

Mechanical stability in such a packing requires that 
there exists no floppy modes, i.e. no collective motions 
of the degrees of freedom of the system (that include the 
Nd degrees of freedom of the particles and changes in 
the box size) for which the distances between objects in 
contact (including both particles and the box) are fixed. 
If such a floppy mode existed, the system would flow 
along it. Rapidly compressed, or poly-disperse packings 
of hard frictionless spherical particles are in fact isostatic 
[Ml HOI li5P5] : the average number of contacts be- 
tween particles, known as the coordination number z, 
is just sufficient to guarantee mechanical stability and 
to avoid the presence of floppy modes, corresponding to 
z = Zc = 2(i [i5H15] . In an isostatic system the removal 





FIG. 2: Illustration of the perturbation considered in this 
work. We select the contact a (left panel), and displace the 
pair of particles that form the contact a by a distance s. 

of any contact leads to the creation of one floppy mode, 
whose properties govern the magnitude of the force found 
in the same contact before its removal, as we shall exem- 
plify below. 

Floppy modes can be generated as follows [44]: two 
particles 1 and 2, forming a contact labelled a, are pushed 
apart while all the other contacts remain closed, as rep- 
resented in Figs, dandjaj We denote by ^^-"^ (s) the dis- 
placement of particle i following the opening of the con- 
tact a by a distance s. This displacement field is uniquely 
defined, because only one floppy mode appears when a 
contact is broken, and exists for sufficiently small open- 
ings s, so as to ensure that no new contacts are formed 
in the system. It satisfies the following equation, that 
embodies the fact that contacts (ij) other than a remain 
closed: 



5R!g>{s)-A,,+0{s^) 



(1) 



where 5, 



1 if and only if the pair (ij) is equal to the 
pair a, and is zero otherwise, and n^j = Rij /rij is the 
director pointing from particle i to particle j. 

On the other hand, force balance in the unperturbed 



packing can be written as: 



Vi, Fi+^f^jrHj 

jii) 



(2) 



where the sum is on all particles j{i) in contact with i, 
Fi is the force exerted by the wall on particle i (and is 
thus zero for particles in the bulk), and fij > is the 
magnitude of the (purely repulsive) force in the contact 
(ij). Multiplying Eq. ^ by any displacement field SRi 
and summing on all particles leads to the virtual work 
theorem: 



i {ij} 



Tlijfi^ 







(3) 



where denotes the summation over all contacts {ij). 

In our system external forces only stem from the bound- 
aries, and the work associated with the displacement field 
5R is Fi ■ 5R, = -pSn. 

Applying Eq. ^ with the vector field Sr['^\s), we 
obtain [T]: 



pSn 



(«)- 



^ • SR. 



(4) 



Introducing the linear fioppy mode 



(a) _ 



dR, 



(a) 



ds 



s=0 



one obtains by derivation: 



(5) 



(6) 



It is convenient to introduce the vector field V*l"' cor- 
responding to the restriction of the fioppy modes to the 



boundary, i.e. V*^f^ = v}"^ for all particles i on the 



boundary, and V' 



otherwise. Then we may write: 



J2f.-v*\ 



(a) 



|F|| • 



(7) 



where = J2i is the angle made between 

the vector fields Fi and Eq. ^ indicates that the 

amplitude of a contact force is governed by how the floppy 
mode associated to that contact couples to the external 
forces at the boundaries. 

A contact force can thus be small for two reasons. 
First, the amplitude of the displacement of the floppy 
mode are smafi, i.e. \\V<°'^\\ < where ()^ de- 

notes averaging over all contacts. For a typical contact 
it was shown that in an isostatic system v}'^^ ~ 1 [29], 
as illustrated in Fig. [Sja. However in some cases the lo- 
cal organization of the particles around the contact a is 
such that local displacements are weakly coupled to the 
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FIG. 3: Examples of two floppy modes, i.e. displacement field resulting fi'om pushing apart a pair of particles carrying a weak 
contact force, represented as shaded. Panel a) exemplified the case where the displacements of the rest of the particles are of 
the same order of the displacements of the pushed pair (corresponding in our notation to fe ~ 1), and the contact force is small 
because the floppy mode is almost orthogonal to a compression of the box. Panel b) displays a different displacement field 
generated in the same configuration as in panel a) by pushing apart a different pair of particles. This time the pushed pair's 
displacement is significantly larger than the displacements in the rest of the system, i.e. & ^ 1. Such contacts are very weakly 
coupled to external stresses applied on the boundary, i.e. they are mechanically isolated. Panel c) displays a local configuration 
of particles that gives rise to small displacements when opening the vertical contact. Even if fo ~ {/), the force in the vertical 
contact can be small if the angle 9 is small, and displacements resulting from opening that contact will be of order 6 ~ sinS. 



rest of the system, as exemplified in Fig. [3jb where sig- 
nificant motion occurs only near the perturbed contact. 
A local model of particle organization generating such 
a weak coupling with the rest of the system is sketched 
in Fig. [3]c. To quantify this behavior we introduce the 
quantity: 



(8) 



ba characterizes the floppy mode in the far field, which is 
of magnitude ^ ba- Note that contacts satisfying 

ba 1 are mechanically isolated: changes of external 
forces applied on the boundaries have little effect on the 
force amplitude, as implied by Eq. ([t]). 

The second possibility is that the floppy mode is 
nearly orthogonal to the external forces, i.e. cos{6a) ^ 
(cos(6'^))^. This situation is illustrated in Fig. [sja. To 
take these possibilities into account it is convenient to 



introduce 



typ 



\\F\\{\\y*^^m)(){cosi90))p, and the di- 



mensionless quantities: 

Wa = -COsi0a)/{cOs{0p))0. (9) 

With these notations Eq. ([7| can be rewritten as: 



ftyp 



baWa . 



(10) 



As it has been shown [T3] that the distribution of con- 
tact forces in packing decay at least exponentially fast 
at large forces, one thus expects that the distributions 
P(b) and P(W) rapidly decay when their arguments are 
larger than the typical values, i.e. for b > 1 and W > 1. 



At small arguments we assume power-law forms: 

P{b) - 6^ for 5 < 1 
P(W) - W^' for < 1 



(11) 



In what follows we shall make the hypothesis that these 
two variables are independent P{b, W) « P{b)P{W) (see 
numerical validation of this hypothesis in Appendix [A| . 
It is then straightforward to compute the behavior of the 
force distribution at low-forces: 



P{f) = J dbdWP{b)P{W)S{f~bW) - /" 



(12) 



From Eq. ( 12 ) we learn that only the smaller of the two 



exponents 9 and 9' can be extracted from the distribution 
of contact forces P{f). 

We have shown in this Section that weak contact forces 
1 3> /a ~ baWa have two possible origins: either the dis- 
placements that result from opening a contact are small, 
i.e. 60, ^ 1, or the same displacements are weakly cou- 
pled to a compressive strain, i.e. Wa ^ 1 (or both). In 
the following we will show how pushing against contacts 
carrying weak forces may lead to rearrangements of the 
contact network, and how the response to this perturba- 
tion is related to the stability of jammed packings. We 
will argue that the exponents 9 and 9' characterize two 
distinct pathways for structural rearrangements, and re- 
late them to several aspects of the geometry and struc- 
ture of packings. 
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III. NONLINEAR STABILITY ANALYSIS 

In this Section we consider the stabihty of the system 
toward rearrangements of the contact network, in par- 
ticular the possibiHty that the volume can be decreased 
by opening a contact and pushing the system along a 
floppy mode (see Figs. [2] and [s]). We consider again the 
displacements (5i?^"^(s) that emanate from pushing apart 
the contact a. Requiring that all other contacts are main- 
tained in tact, and keeping terms up to second order in s, 
we arrive at the following equation for the displacements 
<5^(")(s): 



where 5R^if ■ nj] 



2r, 



+ 0(s3) = s5„ 



(13) 



indicates the projection of onto the 



space orthogonal to Hij. Eq. (13) is analogous to Eq. ([ij), 
in addition to the inclusion of second order terms in s. We 
next apply Eq. M with the displacement field 6R^°'\s) 
defined by Eq. ( [l3| and rearrange, to obtain: 



(iJ) 



We introduce the dimensionless positive number: 

N{f)syao - N{f)/ao 

^ (15) 
and may then rewrite Eq. ([l4| as: 



(16) 



pSn = sU - s^Ca,N{f)/ao + 0{s^) 



As the contact force /„ is positive, Eq. ( 16 ) implies that 
for a sufficiently small opening s the volume always in- 
creases, as illustrated in Fig. Q. However, for larger 
openings the nonlinear term becomes important, and in 
particular if the opening s satisfies: 



s > s 



cM)N ' 



(17) 



a denser packing can be generated [T]. The distance s 
by which a contact a can be opened is bounded how- 
ever, since the displacements of the particles in the entire 
system must eventually lead to the formation of a new 
contact somewhere in the system. We denote by s^ the 
opening created between the pushed pair a when a new 
contact is formed. If < s* , stability is achieved: no 
denser packing can be generated along the direction of 
the floppy mode considered, as illustrated in Fig. [4] 
We define the stability index 
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FIG. 4: Illustration of Eq. (16 1: the energy change pSO, of 



opening a contact a by a distance s. If the resulting displace- 
ment is blocked at an opening of < s* due to the creation 
of a new contact elsewhere in the system, the packing is stable 
against opening the contact a. On the other hand, if the con- 
tact Q can be opened to distances s > s*, a denser packing 
can be found, and thus the packing would be unstable against 
pushing apart the contact a. 



If < 1, the packing is unstable against pushin g ag ainst 
the contact a, as p6V,'-"'> < according to Eq. (16), im- 
plying that a denser packing can be generated 



Eq. ([18) 

can be used numerically to test for the stability of any 
contact in a packing, as shown below. 



IV. SCALING RELATIONS AND MARGINAL 
STABILITY 

The stability of the system requires that Kq, > 1 for all 
contacts a. To unravel the consequence of this constraint, 
we perform a scaling analysis of Eq. ( 18 ) . Firstly, we shall 
use that ^ Waba according to Eq. (10). 

Secondly, the quantity Cq, defined as a normalized sum 



in Eq. ( 15 ), has two contributions. There is a set of 0(1) 
particles very close to the opened contact that satisfy 
ll^/""*!! ~ 1- Assuming that these particles also partic- 
ipate in contacts carrying forces of order (/) (as exem- 
plified in Fig. [sjc) , they should lead to a contribution of 
order 1 /N to . The contribution from the bulk is of or- 
der b^, since far away from the opened contact 1 1^/"' 1 1 ^ 6 
and the majority of contact forces are of order (/). We 
may thus write '-^ 1/iV -|- 6^ (a notation implying that 
Ca ^ l/N lib <t: 1/VN and Cq P otherwise). 

Thirdly, the maximal displacement that the pair a 
can be pushed apart before the creation of a new con- 
tact occurs is governed by the minimal gap /imin between 
particles which are not in contact. In particular, a con- 
tact must be created when the motion of the particles in 
the bulk, of typical amplitude &qS^, is of order of /imin, 



leading to baS' 



The minimal gap follows: 



gih')dh' ~ 1/N , 



(19) 



faO-O 



(18) 



where g{h) is the distribution of gap between particles. 
Assuming a power-law behavior of the distribution of 
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gaps, i.e. g{h) ~ h /a\ ^, one obtains 
{aoN-^)/ba. 



i/ba 



Putting these results together, Eq. ( 18 ) becomes: 



N 



62 



(20) 



It is useful to examine the stability of contacts in the 



(b-W) space; following Eq. (201, contacts a satisfying 



hlN- 



(21) 



are stable against opening. In Fig. [5] a contour plot of 
the joint probability P(6, W) is displayed, together with 
the stability limit given by Eq. ( |21[ ), for various system 
sizes N . 



P{b, W) 
10-' 




FIG. 5: Stability of contacts against opening: the color map 
corresponds to contours having equal values of P{b, W) « 
P{W)P{b), which are assumed to take power-law forms, see 



Eq. (111. The dashed lines are deduced from Eqs. (201 
and (21 1, and correspond to the boundary separating sta- 
ble contacts a having ita > 1, to unstable contacts having 
Ka < 1, for different system sizes A''. 



Imposing the stability of all contacts requires that the 
probability of having an unstable contact in the system 
is much smaller than one, or equivalently: 

P{K)dK < ^ ^ J dWdbP{W)P{b) < ^ (22) 



<1 



Decomposing the integration domain in two parts b ^ 
and b ^ 1/Vn one finds that the integral has two 
positive contributions, and that inequality (22) together 
with Eq. ( 11 ) is satisfied iff: 



7 > 
7 > 



1 



2 
1 



9' 



(23) 
(24) 



Each of the bounds ( 23 1 and ( 24 1 can be seen as resulting 



Mode 1 corresponds to contacts whose floppy modes 
are extended in the system, i.e. 6 1 as illustrated in 
Fig. [3ja. The minimal contact force amplitude /min.i en- 
countered for such contacts in a system of size N is simply 
/min,i ~ Wmin ~ Af-i/(e'+i). The reduction of volume 
stemming from non-linearities is dominated by the bulk, 
where the displacements are of order of the smallest gap 
^min in the system, and Cq ~ 1. Thus following Eq.(20) 



the stability index of the weakest contact of that type is 
K ^ 7V7/(i-7)-i/(i+e')^ and stability requires 7 > 
These modes strongly couple to changes in the forces at 
the boundary, and are expected to dominate the plastic 
response occurring, for example, when a shear stress is 
applied, see Section |Vl) 

Mode 2 corresponds to contacts whose associated 
floppy modes have small displacements in the far field, 
i.e. 5 <C 1 as illustrated in Fig. [3]b,c. Mode 2 contacts 
turn out to be more numerous that mode 1 contacts at 
low forces, and thus dominate P{f) at small /. Their 
minimal force is thus the minimal force in the packing 
and follows /niin.2 ^ &min ^ 7V^i/(^+i). For thcsc con- 
tacts the non-linear effect leading to a reduction of vol- 
ume is dominated by the local displacements near con- 
tact, i.e. c ~ 1/A'', whereas in the bulk the displacements 



are much smaller and of order h„ 



Thus mode 2 corre- 



sponds to a local buckling event of a few particles. This 
local buckling is, however, stabilized by the creation of 
a new contact in the far field. Eq. (201 for the weak- 
est contact of that sort leads to k ~ iV^T/(i-T)-2/(i+9), 
and stability requires 7 > Mode 2 contacts are 

weakly coupled to external forces, and the frequency at 
which they yield under an increasing applied shear stress 
is much smaller that Mode 1, see Discussion. 

In the following we analyze numerically the structure 
and geometry of jammed packings of hard spheres. We 
test the relations between the exponents 7, 0^ 9' given by 
inequalities ( 23 1 and ( 24 1 , and find that they are approx- 



imately saturated. This supports that the distribution 
of contact forces and the pair distribution function are 
coupled at random close packing. In addition, we ver- 
ify the predicted existence of two distinct modes of re- 
laxation, namely mode 1 and mode 2, and demonstrate 
their marginal stability. 



RESULTS FROM COMPUTER 
EXPERIMENTS 



In this Section we numerically test the predictions 



made in Section IV namely the relations between the 



scaling exponents 7, 9 and 9' given in Eqs. (23 1 and (24 1, 



from the stability of a distinct relaxation process. 



and assess the stability of jammed packings. We pro- 
duced an ensemble of jammed packings of various sys- 
tem sizes, and directly measured all of the structural and 
geometrical quantities discussed in the previous sections. 
For a detailed description of the numerical methods and 
calculations used in this work, see Appendix [B| 
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FIG. 6: Distributions of gaps between particles which are not 
in contact, measured in our ensemble of jammed packings 
for different system sizes. The continuous line represent the 
power law g{h) ~ /i"" ''*. Inset: The same distributions of 
gaps measured in jammed configurations before eliminating 
rattlers. The tails of the distributions are different compared 
to the those of the rattler-free jammed packings, and seem to 
follow g{h) ~ /i""'^ (which is represented by the continuous 
line in the inset). This means that rattler- free packings have 
much less small gaps compared to packings with rattlers, as 
previously known [101 lll| . 



A. Scaling relations 

We begin our analyses with the distribution of gaps 
between particles which are not in contact, g{h), which 
is plotted in Fig. [6] This quantity is computed after 
removing the small fraction of rattlers that do not belong 
to the rigid structure. This is the relevant observable in 
our problem, as explained in Appendix B. We observe 
strong finite size effects in the shape of the distributions 
in the small gaps range: the power law behavior seems to 
break down at some system-size dependent gap which is 
decreasing with increasing the system size. In the range 
in which the distributions agree, they obey a power law. 
The divergence of the distribution at small gaps follows 
g{h) ~ ft.""*" with 7 ss 0.38. We note, however, that 
larger jammed configurations are required to improve our 
estimation of this exponent. 

We next turn to the distribution of (see Eq. ([8| 
for definition) which represents the magnitude of the far- 
field displacements of the fioppy mode emanating from 
pushing apart the contact a, see Appendix [B| for details 
about the numerical calculations. Fig. [7] presents the 
distribution P{b) measured in our jammed packings of 
size TV = 1000. We indeed find that P{b) - b'^ follows a 
power law for small 6, with 9 w 0.17. 

In Figure [8] we plot the distribution of the couplings 
Wa of the flonny modes V" to a compressive strain, de- 
fined in Eq. Here too we find that P{W) ^ W^' 
follows a power law for small W, with 9' « 0.44. 

With numerical measurements for the scaling expo- 
nents 7, 9 and 9' in hand, we now turn to validate 



FIG. 7: Distribution P(b) characterizing the mechanical de- 
coupling of contacts. 




FIG. 8: Distribution P{W) characterizing the angle made 
between a fioppy mode and a compressive strain. 




FIG. 9: Distribution P{f) of contact forces / measured in 
our ensemble of jammed packings of systems with A'^ — 1000 
particles, under a pressure p = 1. We do not observe any 
systematic system-size dependence in the same distributions 
measured for different A'^. Inset: the mean minimal contact 
force in a packing vs. packing size A''. 



relations (23) and (24); for Eq. (23) we find 0.38 



7 > (1 -9)/2 w 0.41, whereas for Eq. ^ we find 
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FIG. 10: Stability of contacts plotted in the (b-W) plane - red squares represent contacts that, when perturbed, lead to a 
denser packing, whereas black circles represent stable contacts. The regions of instability and their system-size dependence is 
consistent with our scaling analysis as displayed in Fig.[5j and clearly indicate the existence of two distinct types of instabilities. 
We note that although mode 2 contacts are more numerous than mode 1 contacts (and therefore dominate the distribution of 
contact forces at low force), the number of unstable mode 2 contacts is smaller than the number of unstable mode 1 contacts, 
as manifested in the different prefactors of the scaling of the fraction of unstable contacts, plotted in the rightmost panel. 



0.38 w 7 > 1/(2 + e') w 0.41. The very slight violations 
of the bounds are within the numerical errors, and our 
measures support that the bounds are in fact saturated. 

In Fig. [9] we plot the distribution of the contact forces 
P(/) in our jammed packings. We find (as shown be- 
fore in [TH] for jammed configurations obtained during 
shear flows) that the distribution behaves at small forces 
as P{f) f/{fY^^ with e « 0.18, a power-law that 
appears to hold nicely on three decades. Another test of 
this exponent appears in the inset of Fig. [9] that displays 
the dependence of the mean minimal force with system 
size N. We observe (/min) ~ jv~°'^^, corresponding to 
6 = 1/0.86— 1 « 0.16, in good agreement with our direct 
fit. The se findings validate our prediction spelled out in 
Eq. ( [12I a ccording to which P{f) ~ /min(e,e')^ p^^_ 
pendix \K\ we demonstrate the independence of the vari- 
ables 6 and W , which is assumed to make the prediction 



of Eq. (121 



Note that P{f) and g{h) depend on the protocol by 
which jammed packings are generated, but apparently 
not on the spatial dimension [HI [46^ . For hard sphere 
packings obtained via a Stillinger algorithm Charbon- 
neau et al. [H] found that P{f) - and gih) ~ 

i.e. 7 « 0.42, and to 61 = 0.28 assuming that 



9 < 9'. This allows us to test Eq. 0.42 « 7 > 

(1 — 9)/2 0.36; i.e. a satisfied bound but again quite 
close to saturation. 

For soft decompressed particles these authors found 
P{f) ^ /°'*^ and 7 w 0.39. E 9 < 9' the bound of 
Eq. ([23| gives 0.39 w 7 > (l-6')/2 w 0.29, i.e. a satisfied 
bound but now somewhat further away from marginality. 
Note however that since the exponent they find for the 
force distribution is large in this case, and close to our 
value of 9' (see below), we suspect that the condition 
9 < 9' does not hold in this case. Ii9' < 9 then the bound 
Eq. (24) should be used instead, and reads 0.39 « 7 > 
1/(2 -I- 9') 0.41, i.e. consistent with the marginality of 
mode one, as we also found with very similar exponents. 



B. Distinct modes of instability 

With the measurements of sj^, Cq., and Wa in hand, 
we can directly measure the stability index of each 
contact a in each jammed packing of our ensemble via 
Eq. (18 1. The stability of contacts on the (b-W) plane 
is visualized in Fig. [lO[ We find a consistent picture to 
that described in Fig. [51 Two distinct modes of relaxation 
emerge from this analysis: mode 1, which are contacts a 
having 6q, ~ 1 and Wa <C 1, and mode 2, which are 
contacts a having ba <^ 1 and Wa ~ 1. This result 
demonstrates that contacts carrying small forces /„ ^ 
baWa are potentially unstable, though the mechanical 
nature of the instability can belong to either one of the 
two distinct modes. 

We found that roughly two contacts in a packing, if 
opened, may lead to a denser packing, and are therefore 
unstable. Accordingly, we find that the fraction of unsta- 
ble contacts, plotted in the rightmost panel of Fig. (10 1, 



scales as 1/N for Mode I and II respectively, indicat- 
ing that our packings are indeed marginally stable. This 
further vindicates our findings that the inequalities pre- 
sented in Eqs. (23) and (24) appear to be saturated. 



VI. AVALANCHES OF PLASTIC EVENTS 



When Eqs. (23) and (24) are strictly satisfied, for 



large packings one never observes instabilities: in that 
case <C s* even for the contacts carrying the weakest 
forces, implying according to Fig. [4| that opening a gap 
between particles simply leads to a linear increase of vol- 
ume. However, empirically we observe that the bounds 
of Eqs. (23) and (24) appear saturated, and as shown in 



Fig.[l0|we find that there are instabilities in our packings, 
both for type 1 and type 2 contacts. Thus if temperature 
is present, or if shear stress is applied, such instabilities 
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will occur, and one must consider the rate at which these 
events take place, and their consequences. 

We shall consider the case where external forces (corre- 
sponding, e.g., to a shear stress) are applied. According 
to Eq. (l6| applying an external force AFi on some par- 
ticle i changes the contact force in some contact a by an 
amount 1^: 



A/„ = ■ AF, ^ b^AF, 



(25) 



For a contact to yield its force must vanish, i.e. A/q, < 
—fa — baWaftyp, thus the threshold stress where a rear- 
rangement occurs is AF ~ Waftyp- it is large for type 2 
contacts for which Wa ^ 1, but very small for type 1 con- 
tacts where Wa can be vanishingly small. This argument 
supports that mode 1 rearrangements are much more fre- 
quent and presumably govern the plastic response under 
an imposed stress. 

Let us assume that mode 1 and/or mode 2 are 
marginally stable, and consider the consequence of their 
respective relaxation. Let /3 label the contact that closes 
at . By symmetry, the results obtained when the con- 
tact a was opened and /? was closed also apply to the 
newly obtained configuration if /3 is re-opened and a is 
re-closed. In particular, the relation 



becomes 



_d_ 



dsi 



[pSV{sa)] 



(26) 



(27) 



where is the distance by which the contact /3 is opened. 
Since both for modes 1 and 2, the contact /3 is far from the 
contact a, one has dsp/ds^ ~ &a- Furthermore, marginal 



stability corresponds to s' 
This implies that 



dSa 



[p5V{s^)] 



from which we get: 



dSa 



for the weakest contacts. 



[pSV{sa)] , (28) 



d 

dSa 



[pSVisa)] 



(29) 



Thus the force in the new contact is of order of the con- 
tact that was opened for mode 1 relaxation for which 
b ~ C(l), but it is much larger for mode 2, for which b 
can be vanishingly small. Forming a new contact is equiv- 
alent to applying a dipole of external force of the system. 
According to Eq. (251, this will change the contact force 
in any contact 6 by some amount Afs ^ faibs/ba). 

Consider that a type 1 instability occurs, i.e. /„ ^ 
/i.min and ba ^ I- The effect on the weakest type 2 
contacts, labeled 5 = (2, min) is negligible as in that case 



A/2 



fi 



h 



ftyp 



-fs < /2 



(30) 



However the weakest type 1 contacts has a finite probabil- 
ity to become unstable and to open, as A/i^min ~ /i,min- 
Thus when the system is marginally stable, the relax- 
ation of a Mode 1 contact can lead to a sequence of 
events where other such contacts become unstable and 
open. Numerically such avalanches of relaxation are seen 
and are power-law distributed f415, a situation that is 
only possible if marginal stability is satisfied [1,. 

The same argument applied when a type 2 contact re- 
laxes leads to a similar conclusion: the marginal stability 
of Mode 2 corresponds precisely to the situation where 
the relaxation of a type 2 contact can trigger the insta- 
bility of another such contact with a finite probability. 
However one finds a key difference: if a mode 2 contact 
relaxes it will affect significantly all the contacts forces 
in the system, thus resulting in a large global change of 
the contact network. Thus Mode 2 relaxation appears to 
be extremely rare and hard to trigger, but to have a very 
dramatic effect when it occurs. 



VII. SOFT PARTICLES 

Packings of hard particles are isostatic, and the re- 
sponse to a local strain does not decay with distance, 
unlike generic amorphous materials in which elastic re- 
sponses to local strains decay as power laws away from 
the perturbation. It is thus important to consider how 
our results extend to more realistic models, for example 
to soft compressed particles (</> > 0c) or thermal colloidal 
systems (0 < 0c)- Do the relaxations processes that we 
have found in packings of hard particles still exist in such 
systems, and if so what fixes their density? 

Previous works [71 [51 [311 [IZ] aiming at relating relax- 
ation processes to the microscopic structure in these sys- 
tems have focused on the linear response, in particular 
on the soft vibrational modes very different from plane 
waves that are present at low-frequency. The presence 
of such modes is expected if these systems live close to 
an elastic instability, as is indeed observed after numeri- 
cal quenches into the glass phase are performed, at least 
for a rather large interval of packing fractions around 0c 
O [28l HQ] . The idea that nearly unstable modes would 
generate low- activation barriers above which the system 
can relax is rather natural, and is supported by numerical 
observations showing that thermal relaxation mostly oc- 
curs along a few soft modes both in colloidal systems [7] 
and in Lennard- Jones systems [35] . An accurate descrip- 
tion of the barriers associated with soft modes is however 
lacking j49j and would be very interesting to investigate 
numerically. 

The present work suggests that when the interaction 
potential is non-linear (for example soft spheres where 
the potential is a one-sided harmonic spring), intrinsi- 
cally non-linear relaxation processes should be consid- 
ered too. This view that the relaxation processes of 
type 1 and/or type 2 play a role away from the jamming 
threshold is supported by the recent numerical results 
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of Charbonneau et al. [T2], where P{f) and g{h) have 
been carefully measured as the packing fraction is var- 
ied, both in colloidal systems {(f) < (pc) and in systems 
of soft particles > 0c)- For concreteness, we shall fo- 
cus on soft particles, which in that study are harmonic, 
and where the stiffness of the interactions is taken to be 
unity. These authors find that under a compression of 
amplitude A(/>, power-law tails are still observed in the 
distribution functions of structural quantities, which in 
our notation writes 



PU) 



if) J 



and g(h) 



ao 



(31) 

excepted on an extremely narrow interval of gaps h < 
h* ~ Acpt" and forces f < f* ^ with /i w 1.75, 

where these distributions do not significantly vary. As 
pointed out in [12', the exponent fi can be related to 7 
and to the force exponent (imii{9,0') in our notation). 
To see this, let us generalize the definition of the gap 
distribution function g{h) to include for the possibility 
of negative gaps, i.e. particles in contact for which = 
— ha- Assuming that g{h) ~ l/h"^ for h ^ A(f>^ and 
g{h) - (-/i)«""(«.»')/A0™"(».«') (as required by the force 
distribution) for —h ^ —Acf)'^, and that g{h) smoothly 
varies for \h\ <^ Acj)'^ implies by continuity that /x = (1 + 
mm{e,e'))/{j + inm{e,e')). li{9 < 9') marginal stability 
of mode two (see Eq. (23)) implies ^ = 2, whereas if 
9' < 9 marginal stability of mode one would imply /i = 
(1/7 — l)/(7 -I- 1/7 — 2) « 1.64, in closer agreement with 
the observations of [12) . 

The observations that force and pair distributions arc 
nearly unchanged in compressed packings, excepted on a 
very narrow range of forces and gaps, suggest that the 
modes we introduced still play a role in controlling the 
structure in this situation. The numerical data of [12 
are consistent with the idea that occurrences of mode 2 
are rare in compressed packing, perhaps because they are 
too fragile and "buckle" under compression. To test this 
scenario and distinguish which of the exponents 9 or 9' 
is smaller in soft spheres, one would need to perform a 
numerical analysis similar to the one presented in this 
work. 



VIII. CONCLUSION 

We conclude by a summary of our results and open 
questions. We have established the existence of two dis- 
tinct non-linear modes of relaxation in random packings 
of hard frictionless spheres, associated with changes of 
the network of contacts. Mode two corresponds to the 
local buckling of a few particles. Such excitations are 
numerous and their density is described by an exponent 
9 that also characterizes the distribution of contact forces 
P{f) ^ ■ Using a stability argument we proved that 
the density of such modes is bounded from above, lead- 
ing to 7 > -'-f^, where 7 characterizes the singularity 



of the pair distribution function g{r) at contact. Mode 
one, whose density is characterized by the exponent 9' , 
corresponds to extended motion of the particles, and the 
bound on their density is 7 > 2+¥' ^® performed nu- 
merics that support that these bounds are saturated with 
7 = 0.38, 9 = 0.17 and 9' ^ 0.44. In addition, we intro- 
duced numerical methods to analyze the entire sets of ex- 
citations of a given packing, which indicate that the two 
modes of relaxation considered are indeed barely stable. 
These results support that marginality is the fundamen- 
tal principle governing which ensemble of configurations 
are visited in this out-of-equilibrium system, and con- 
strain key aspects of the microscopic structure of pack- 
ings, such as the distribution of force P{f) at low forces 
and the pair distribution function g(r) at small gaps. It 
also yields a natural explanation for the crackling dy- 
namics observed in driven packings of hard frictionless 
particles. 

Concerning the random close packing of hard particles, 
three questions stand out. (i) If marginal stability is 
indeed found, as appears to be the case for our system 
preparation, the three exponents 9, 9' and 7 are related 
by two relationships. What then, fixes the exponent 7? 
(m) What explains the difference between soft and hard 
spheres, and are mode two marginal or even present in 
soft spheres as well? (Hi) Hard sphere packings display 
power-law distributed avalanches of plasticity 41j. Can 
the exponent characterizing the size of the power-law be 
computed? 

to ex- 



Next it is important, as discussed in Section VII 



tend this approach to situations where particles are elas- 
tic and compressed or to colloidal glasses. The structure 
still appears to be sensitive to the relaxations modes we 
have introduced. Is there still some kind of marginality 
present? What fixes the concentration of these relaxation 
processes? 

Lastly, there is evidence that these relaxation processes 
play an important role in flow. Anisotropic configura- 
tions that jammed during flow display a distribution of 
force nearly identical to isotropic packings [13'. Away 
from the jamming threshold this distribution is cut off 
below some relative force scale that vanishes near jam- 
ming [131, similarly to what happens in soft particles or 
colloidal systems [H] . These results suggest to develop a 
description of contact dynamics and its effect on the rhe- 
ological properties of dense suspensions based on these 
excitations l50l. 
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Appendix A: Independence of b and W 

The quantities b and W, defined in Eqs. (|8| and ([9| 
respectively, are assumed to be statistically independent 
in the range b 1 and W 1. This assumption is 
used when the prediction of the force distribution form is 
made in Eq. ( |12[ ), and in the scaling arguments made for 
Eq. ( 22 ) . In this Appendix we demonstrate how well this 



assumption holds by analyzing data from our numerical 
simulations. For a detailed description of the numerical 
methods and of calculations of structural quantities, see 
Appendix [B] 




FIG. 11: Conditional distributions: a) P(VF|fe), b) P{b\W). 
Both conditional distributions seem to converge to a limit- 
ing distribution when the conditioned variable is smaller than 
unity. 

If b and W were independent, one would observe that 
the conditional probabilities P{b\W) and P(VF|6) obey 



P{b\W) = P{b) and P{W\b) = P{W) 



(Al) 



namely that the conditional distribution does not depend 
on the value of the conditioned variables. In Fig. [TT] we 
exemplify that for & << 1 and 14^ << 1 this is indeed the 
case. We conclude from these results that in these limits, 
the joint distribution P{b, W) « P{b)P{W) - b'^W'^' , see 
Eq. ITTl). 



Appendix B: Numerical methods 

To generate numerically jammed packings under pe- 
riodic boundary condition we used the event-driven 



method we developed in [5T]; the method enables one 
to isotropically compress systems of hard frictionless 
spheres under overdamped dynamics. We initialized each 
independent compression run with a random initial con- 
dition at a packing fraction (/iinitiai < 4>c ~ 0.645 lower 
than the critical packing fraction 0c, and compressed 
the systems until they jammed. We generated about 
1000 jammed configurations for N = 2000, about 4500 
jammed configurations for N — 1000 and about 9000 
jammed configurations for N = 125. In the compression 
simulations wc bound the mean error in contacts over- 
lap to be smaller than 10~®ao |51| . where uq is the mean 
diameter of the particles. We employed a 50:50 binary 
mixture of large and small particles, such that the ratio 
between the radii of the large and small particles is 1.4. 

In jammed configurations, a small percentage of the 
particles are loosely connected to the rigid backbone of 
the packing [51]. We eliminated these 'rattler' particles 
by fixing the position of the all the particles belonging to 
the rigid backbone, and applying a 'gravitational' force 
field such that the rattlers sink to the bottom of their re- 
spective cavities, and create at least three contacts with 
neighboring particles. We generalized the same method 
described in [21] for this step of our numerics. Eliminat- 
ing the rattlers turns out to be important to study quanti- 
tatively the stability of packings - although how they are 
removed, by sinking or by taking them out of the system, 
leads to the same results. Packings in which rattlers are 
not eliminated posses many more small gaps compared 
to rattler-free packings. These small gaps would affect 
our numerical estimation of the stabilization distances 
where new contacts are formed. Thus, the amplitude 
of the displacements SR^"^ resulting from pushing apart 
contacts would appear more limited than what they ac- 
tually are, since new contacts created with rattlers do 
not confer mechanical stability. 

The structural quantities of interest for the analysis 
described in this work are the contact forces /„ defined 
in Eq. the magnitude of the far-field displacements 
bry defined in Eq. ([sl, the cosine of the angles Wa as 



defined in Eq. ([9|, and the stabilization distances sj, at 
which a new contact is formed. The extraction of these 
quantities from our jammed packings is explained in the 
following buUcted paragraphs: 

• Contact forces are resolved by solving Eq. ([2|, to- 
gether with the constraint that the pressure is fixed: 



— pfld 



(Bl) 



• Calculation of the far-field displacements ba re- 
quires the availability of the floppy mode V^"' 
which emanates from pushing against a contact a, 
defined in Eq. A jammed packing embedded 
in a fixed container possesses no such modes, and 
in the arguments of Sect. |IV| we let the walls con- 
fining the system move, thus allowing for the vol- 
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ume to vary. The same situation occurs with pe- 
riodic boundary conditions, where the system size 
needs to be aUowed to change length to generate 
one floppy mode by opening one single contact. We 
thus let the periodic system change length, or met- 
ric, and define e(")(s) = (L(s) - L(0))/L(0), the 
compressive strain that occurs as the contact a is 
opened by s. In this geometry the floppy mode 
equation is not Eq. (13), but rather: 



2r, 



-f e(")(s)r,, 



(B2) 

From Eq. (B2) it is straightforward to derive our 



previous results that remain unchanged, in partic- 
ular Eqs. (14) and (18). Numerically, we solve the 
equation 



(B3) 



for the floppy mode V^"> associated with the open- 
ing of any contact a. With the floppy mode V^°'^ in 
hand, we can directly calculate the dimcnsionless 
constants Ca as defined in Eq. (15). 



The definition of ba proposed in Eq. ^ is not 
practical in a periodic system, and numerically we 
choose the following definition: 



6a = median f v}"''^ \ , (B4) 

contacts (ij) L J 



where the median is taken over all contacts. It is 
important to consider the median and not the mean 
of the relative displacements as the mean is affected 
by the few particles near contact for which V ^ 1, 
and thus does not represent the amplitude of dis- 



placements in the far field. Fig. 12 shows the distri- 
bution of relative displacements following the open- 
ing of a contact, and exemplifies its dependence on 
the strength of the contact force. These observa- 
tions demonstrate that when a contact a is opened 
and displacement along a fioppy mode occurs, a 
well-defined displacement scale ba emerges. 



a) 10- 




10.000 10.000 
contact # 



FIG. 12: a) Scatter plot of relative displacements of particles 
in contact, following the pushing apart of a contact carrying 
a weak force fa ~ 10~^ (crosses) and a typical force /„ « 1 
(dots), both measured in a single jammed configuration in two 
dimensions with TV = 10000 particles. Notice that even in the 
weak force case (crosses), there are a few contacts with rela- 
tive displacement of order unity, which is why considering the 
median of relative displacements (and not the mean) is use- 
ful, b) Distributions of relative displacements rescaled by the 
medians be, defined in Eq. (B4l. The symbols are consistent 
with the left panel. 



• Once the far-field displacements ba are calculated, 
then following Eq. ( [lO| the cosine of the angles Wa 
of Eq. (|9| are calculated as 



Wa 



fa/ba ■ 



(B5) 



The stabilization displacements s\. are the maxi- 
mal distance s that the pair a can be displaced 
before creating a new contact. To calculate sj,, we 
consider all pairs of particles {ij) that are not in 
contact, and use the scheme for finding the next 
collision time in a hard sphere gas |52j . treating 
V^"^ as velocities and s as time. 
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